Índice
Nova shape¶
Imports¶
In [1]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as cls
from pandas import Series, DataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes, mark_inset
import matplotlib.ticker as ticker
from collections import defaultdict, namedtuple
from pandas import ExcelWriter
import xlsxwriter
import pandas as pd
import calendar
import pickle
import datetime
from scipy import stats
import itertools
import statistics
# set seaborn plot style
import seaborn as sns
sns.set()
sns.set_context("paper")
sns.set_style("whitegrid")
# set default font
plt.rcParams["font.family"] = "Times New Roman"
In [7]:
def get_ecos():
m = Basemap(projection='cyl')
shpf = m.readshapefile('C:/Users/alpha/Documents/MODIS/MCD64A1/newreprojected/newreprojected','ecos',linewidth=0.1)
return m.ecos_info, m.ecos
In [6]:
m = Basemap(projection='cyl')
shpf = m.readshapefile('C:/Users/alpha/Documents/MODIS/MCD64A1/newreprojected/newreprojected','ecos',linewidth=0.1)
ecos_info, ecos = m.ecos_info, m.ecos
In [4]:
%qtconsole
In [10]:
ecos_info[0].keys()
Out[10]:
In [11]:
months = [calendar.month_name[i +1][:3] for i in range(12)]
In [12]:
years = [str(y) for y in range(2001,2018)]
In [42]:
ecoids = [e['ECO_ID'] for e in ecos_info]
bionums = [e['BIOME_NUM'] for e in ecos_info]
bionames = [e['BIOME_NAME'] for e in ecos_info]
ids = [e['ID'] for e in ecos_info]
climes = [e['GRIDCODE'] for e in ecos_info]
continents = [e['CONTINENT'] for e in ecos_info]
countries = [e['FIPS_CNTRY'] for e in ecos_info]
areas= [e['Area_polig'] for e in ecos_info]
In [59]:
for i ,id in enumerate(ecoids):
if id == 0:
bionums[i] = 0.0
In [20]:
result = {}
for i in range(1,18):
suf = '20{0:02d}'.format(i)
in_file= 'newresult' + suf
result[suf] = DataFrame(np.load(in_file), columns =months)
In [60]:
for i in range(1,18):
suf = '20{0:02d}'.format(i)
result[suf] = result[suf].assign(eco_id=ecoids,biome_num =bionums,biome_name=bionames,idd=ids,clime=climes,continent=continents,country=countries, area = areas)
In [92]:
grouped = result['2001'][ ['eco_id','continent','clime','country','idd' ,'area'] ].groupby(['eco_id','continent','clime','country','idd' ,'area'])
areas_row = []
for a, _ in grouped:
areas_row.append(a[-1])
In [105]:
attrs_table = {}
for yr, df in result.items():
attrs_table[yr] = df[ months + ['eco_id','continent','clime','country','idd','biome_num' ,'area']].groupby(['eco_id','continent','clime','country','idd','biome_num','area']).sum()
In [110]:
at_table = pd.concat(attrs_table)
Excel¶
In [111]:
writer = ExcelWriter('newshape.xlsx',engine='xlsxwriter')
In [112]:
at_table.to_excel(writer,sheet_name='Burned Pixels',engine='xlsxwriter')
In [113]:
writer.save()
In [206]:
plt.rcParams['mathtext.fontset'] = 'custom'
plt.rcParams['mathtext.rm'] = "Times New Roman"
plt.rcParams['mathtext.it'] = "Times New Roman"
plt.rcParams['mathtext.bf'] = "Times New Roman"
In [207]:
g = sns.regplot(np.arange(12),attrs_table['2002'].sum()/1e6, fit_reg= False)
ylabels = g.get_yticklabels()
g.set_xticklabels(months[::2], fontsize=15)
g.set_xticks(np.arange(0,12,2))
g.set_ylabel('Burned Pixels $(10^6)$', fontsize=15)
plt.setp(ylabels,fontsize=15)
#g.set(xticklabels=months[::2],xticks=np.arange(0,12,2), ylabel='Burned Pixels $(10^6)$', ylabelsize=15)
plt.show()
In [210]:
plt.plot(np.arange(12),attrs_table['2002'].sum()/1e6)
g =plt.gca()
ylabels = g.get_yticklabels()
g.set_xticklabels(months[::2], fontsize=15)
g.set_xticks(np.arange(0,12,2))
g.set_ylabel('Burned Pixels $(10^6)$', fontsize=15)
plt.setp(ylabels,fontsize=15)
#g.set(xticklabels=months[::2],xticks=np.arange(0,12,2), ylabel='Burned Pixels $(10^6)$', ylabelsize=15)
plt.show()
In [211]:
plt.scatter(np.arange(12),attrs_table['2002'].sum()/1e6)
g =plt.gca()
ylabels = g.get_yticklabels()
g.set_xticklabels(months[::2], fontsize=15)
g.set_xticks(np.arange(0,12,2))
g.set_ylabel('Burned Pixels $(10^6)$', fontsize=15)
plt.setp(ylabels,fontsize=15)
#g.set(xticklabels=months[::2],xticks=np.arange(0,12,2), ylabel='Burned Pixels $(10^6)$', ylabelsize=15)
plt.show()
In [215]:
at_table.shape
Out[215]:
In [221]:
at_table.loc['2001',0,'Antarctica']
Out[221]:
In [222]:
at_table.sum()
Out[222]:
Variação intra-anual. Totais¶
In [246]:
fig, ax = plt.subplots()
ax.plot( np.arange(12),at_table.sum()/1e6)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months, fontsize=15)
ax.set_title('Intra-annual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels $(10^6)$", fontsize = 15)
plt.tight_layout()
plt.show()
Variação inter-anual. Totais¶
In [239]:
at_table.index.names = ['year', 'eco_id', 'continent', 'clime', 'country', 'idd', 'biome_num', 'area']
In [244]:
at_table.groupby('year').sum().sum(axis=1)
Out[244]:
In [245]:
fig, ax = plt.subplots()
labels = years[::2]
ax.plot(years, at_table.groupby('year').sum().sum(axis=1)/1e6)
ax.set_xticks(labels)
ax.set_xticklabels(labels, fontsize=15)
ax.set_title('Interannual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels ($10^6$)", fontsize = 15)
plt.tight_layout()
plt.show()
Variação intermensal. Totais¶
In [249]:
at_table.groupby('year').sum().values.flatten()
Out[249]:
In [251]:
ym = itertools.product(years,months)
In [262]:
ymlabels = [ m + ' ' + y[-2:] for y, m in ym]
fig, ax = plt.subplots(figsize=(16,12))
labels = ymlabels[::6]
ax.plot( np.arange(12*17),at_table.groupby('year').sum().values.flatten()/1e6)
ax.axhline(y=statistics.mean(at_table.groupby('year').sum().values.flatten()/1e6), color='r', linestyle='--', label='mean')
ax.set_xticks(np.arange(0,12*17,6))
ax.set_xticklabels(labels,fontsize=12, rotation=30)
labels_y = ax.get_yticklabels()
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels ($10^6$)", fontsize = 15)
ax.set_title('Intra-monthly Variability',fontsize=15)
plt.legend(fontsize=15)
#prop={'size': 6}
plt.show()
Variação intra_anual em cada bioma¶
In [264]:
at_table.groupby('biome_num').sum()
Out[264]:
In [283]:
at_table.loc['2001',0, 'Antarctica',61,'AY',2250 ]
Out[283]:
In [291]:
areas = at_table.loc['2001'].index.get_level_values(6)
In [296]:
at_table.index.get_level_values(7)
Out[296]:
In [297]:
141338/17
Out[297]:
In [303]:
all_areas =at_table.index.get_level_values(7)
In [312]:
at_table_cp = at_table.copy()
In [304]:
#at_table_cp.areas = all_areas
In [313]:
at_table_cp.reset_index(inplace=True)
In [325]:
agg_dict = {month : 'sum' for month in months}
In [326]:
div_17 = lambda x: sum(x)/17
In [330]:
agg_dict['area'] = div_17
In [333]:
biomes_table = at_table_cp.groupby('biome_num').agg(agg_dict)
In [341]:
biomes_table[months].div(biomes_table['area'], axis =0)
Out[341]:
In [338]:
biomes_table['area']
Out[338]:
In [345]:
# for serialization:
tabelas = {'multindex': at_table, 'simpleindex':at_table_cp}
In [346]:
with open('tabelas','wb') as f:
f.write(pickle.dumps(tabelas))
In [369]:
biomes_intra = biomes_table[months].div(biomes_table['area'], axis =0)
In [356]:
biome_names_dict = dict( (num, name) for num, name in zip(bionums, bionames))
Mesma escala¶
In [372]:
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = months
ax = ax.ravel()
for i in range(15):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, biomes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
#ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
Sem escala comum¶
In [373]:
fig, ax = plt.subplots(5,3, figsize =(18,12))
labels = months
ax = ax.ravel()
for i in range(15):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, biomes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
#ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
Variação interanual em cada bioma¶
In [393]:
at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months]
Out[393]:
In [390]:
biomes_table['area'].size
Out[390]:
In [397]:
## test 2001 may bionum = 5 assert 3.547547e-04
bp = at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().iloc[5]['May']
area_5 = biomes_table['area'][5]
assert bp/area_5 == 3.547547e-04, print((bp/area_5 - 3.547547e-04 )*100/3.547547e-04)
In [399]:
at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months].sum(axis=1)
Out[399]:
In [400]:
biomes_inter = [at_table_cp[at_table_cp.year == yr].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months].sum(axis=1) for yr in years]
In [405]:
biomes_inter = DataFrame(biomes_inter).T
In [406]:
biomes_inter
Out[406]:
Mesma escala¶
In [414]:
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
ax = ax.ravel()
for i in range(15):
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), biomes_inter.iloc[i])
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, biomes_inter.iloc[i], 'o', label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
Sem escala comum¶
In [416]:
fig, ax = plt.subplots(5,3, figsize =(18,12) )
labels = years[::2]
ax = ax.ravel()
for i in range(15):
#slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), biomes_inter.iloc[i])
#ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, biomes_inter.iloc[i], label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(biome_names_dict[i],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout()
plt.show()
In [417]:
biomes_inter
Out[417]:
In [419]:
at_table_cp[at_table_cp.year == '2001'].groupby('biome_num').sum().div(biomes_table['area'], axis=0)[months].sum()
Out[419]:
Variação intra-anual em cada clima¶
In [422]:
no_climas = len(set(climes));no_climas
Out[422]:
In [427]:
climes_table = at_table_cp.groupby('clime').agg(agg_dict)
In [428]:
climes_table[months].div(climes_table['area'], axis =0)
Out[428]:
In [430]:
climes_intra = climes_table[months].div(climes_table['area'], axis =0)
Mesma escala¶
In [431]:
climanames = [e['Clima_Name'] for e in ecos_info]
In [434]:
climes_intra.index[0]
Out[434]:
In [432]:
clima_names_dict = dict( (num, name) for num, name in zip(climes, climanames))
In [466]:
# 11 x 3, 31 climas
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) , sharey=True)
labels = months
ax = ax.ravel()
for i in range(31):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, climes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climes_intra.index[i]],fontsize=15)
#ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
Sem escala comum¶
In [467]:
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) )
labels = months
ax = ax.ravel()
for i in range(31):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(months, climes_intra.iloc[i], 'o')
ax[i].set_xticks(np.arange(12))
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climes_intra.index[i]],fontsize=15)
#ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
In [471]:
clima_names_dict.values()
Out[471]:
Variação inter-anual em cada clima¶
In [488]:
at_table_cp[at_table_cp['year']=='2001'].groupby('clime').sum().div(climes_table['area'], axis=0)[months].sum(axis=1)
Out[488]:
In [493]:
climas_inter = [at_table_cp[at_table_cp['year']==yr].groupby('clime').sum().div(climes_table['area'], axis=0)[months].sum(axis=1) for yr in years]
In [495]:
climas_inter = DataFrame(climas_inter).T
In [496]:
climas_inter
Out[496]:
Mesma escala¶
In [500]:
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) , sharey=True)
labels = years[::2]
ax = ax.ravel()
for i in range(31):
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), climas_inter.iloc[i])
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, climas_inter.iloc[i], 'o', label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climas_inter.index[i]],fontsize=15)
ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
Escala comum¶
In [501]:
fig, ax = plt.subplots(nrows=11,ncols=3, figsize =(18,40) )
labels = years[::2]
ax = ax.ravel()
for i in range(31):
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), climas_inter.iloc[i])
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, climas_inter.iloc[i], 'o', label='')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(clima_names_dict[climas_inter.index[i]],fontsize=15)
ax[i].legend(fontsize=15)
for i in range(31,33):
ax[i].set_visible(False)
fig.tight_layout()
plt.show()
Variação inter-anual por mês em cada bioma¶
In [530]:
groupedby_year_clima = at_table_cp.groupby(['clime','year']).sum()
In [531]:
groupedby_year_clima = groupedby_year_clima[months].div(groupedby_year_clima.area, axis =0)
In [537]:
groupedby_year_clima =groupedby_year_clima.reset_index()
In [544]:
groupedby_year_clima[groupedby_year_clima['clime'] == 11][months].T.iloc[0]
Out[544]:
In [567]:
labels = years[::2]
for climcode in clima_names_dict.keys():
df = groupedby_year_clima[groupedby_year_clima['clime'] == climcode][months].T
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
ax = ax.ravel()
fig.suptitle(clima_names_dict[climcode],fontsize=15)
for i in range(12):
data = df.iloc[i]
slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), data)
ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, data, 'o', label='_nolegend_')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(calendar.month_name[i +1],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
In [569]:
labels = years[::2]
for climcode in clima_names_dict.keys():
df = groupedby_year_clima[groupedby_year_clima['clime'] == climcode][months].T
fig, ax = plt.subplots(4,3, figsize =(20,12) )
ax = ax.ravel()
fig.suptitle(clima_names_dict[climcode],fontsize=15)
for i in range(12):
data = df.iloc[i]
#slope, intercept, r_value, p_value, std_err = stats.linregress(range(17), data)
#ax[i].plot(years, intercept + slope*range(17), 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i].plot(years, data, 'o', label='_nolegend_')
ax[i].set_xticks(labels)
ax[i].set_xticklabels(labels,fontsize=15)
labels_y = ax[i].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i].set_title(calendar.month_name[i +1],fontsize=15)
ax[i].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
In [592]:
# as linhas seguintes dão o mesmo resultado
at_table_cp.groupby(['year','clime']).sum()[months];
pd.pivot_table(at_table_cp, values=months, index=['year', 'clime'],aggfunc=np.sum)[months];
In [593]:
notincluded13 = pickle.loads(open('notincluded2013','rb').read())
df_notinc13 = DataFrame(notincluded13)
df_notinc13.to_csv('not_included13.csv', header=False)
In [610]:
notincluded17 = pickle.loads(open('notincluded2017','rb').read())
df_notinc17 = DataFrame(notincluded17)
df_notinc17.to_csv('not_included17.csv', header=False)